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Introduction 


This  project  was  originally  a  two-year  program  to  develop  unique  microfluidic  modeling 
software  with  a  specific  focus  on  the  suspended  species  in  solution  and  the  unique  near-field 
interactions  that  occur  in  microenvironments,  more  specifically  in  microfluidic  systems.  In 
addition  this  effort  focused  on  the  inclusion  of  external  fields  and  in  particular  AC-  electrokinetic 
effects.  A  flow  diagram  or  road  map  describing  the  major  goals  of  the  effort  is  given  below  in 
Figure  1. 
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Figure  1.  The  original  road  map  for  proposed  effort.  The  light  blue  boxes  represented 
near  term  milestones  and  the  dark  blue  boxes  represented  longer  term  milestones. 


This  road  map  was  more  than  merely  deseribing  the  proposed  effort;  In  a  greater  sense,  it  was 
proposed  to  and  shown  to  give  a  larger  pieture  to  the  needs  of  modeling  in  mierofluidies.  To 
develop  these  unique  eomputational  eapabilities,  we  developed  3D,  parallel  lattiee  Boltzmann 
simulation  eapabilities  [1]. 

The  Lattice  Boltzmann  (LB)  is  uniquely  suited  for  studying  the  dynamic  behavior  of 
macromolecules  in  microfluidic  devices.  Current  methods  -  Finite  Element  Analysis  (FEA)  and 
Boundary  Element  Methods  (BEM)  -  have  been  useful  for  modeling  and  simulating  pure  fluid 
flows  in  microfluidic  networks  with  coupled  heat  transfer  and  limited  electrokinetic  phenomena. 
The  current  microfluidics  simulation  tools  do  not,  however,  explicitly  take  into  account  the 
behavior  of  the  suspended  phase  (macromolecules).  While  valuable,  these  approaches  are 
limited  in  crucial  ways;  namely,  it  is  very  difficult  to  include  finite  solutes  with  complex 
geometries  in  EEA  methods.  BEM  can  easily  incorporate  mobile  particles  with  external  field 
effects;  however,  this  approach  is  only  valid  for  Newtonian  flows  in  either  the  zero  or  infinite 
Reynolds  number  flow  regimes.  As  a  consequence,  it  is  not  possible  to  study  non-Newtonian 
effects  or  small  but  finite  Reynolds  number,  e.g.,  0  <  Re  <  10,  flow  effects  that  can  occur  in 
microfluidic  devices.  The  lattice  Boltzmann  method  does  not  suffer  from  these  limitations.  The 
strength  of  the  EB  method  is  that  it  is  extremely  flexible.  More  specifically,  it  is  very  easy  to 
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incorporate  complex  boundaries  such  as  curved  bounding  walls,  obstacles  and  mobile  solutes 
with  complex  shapes.  The  only  limitation  is  the  user’s  ability  to  mathematically  describe  the 
geometry.  The  LB  method  easily  permits  the  inclusion  of  new  physics.  The  capability  permits 
the  inclusion  of  new  physics,  such  as,  external  AC  and  DC  field  effects  (dielectrophoresis  and 
electrokinetic  phenomena)  and  colloidal  interactions  (electrostatic  and  van  der  Waal  interactions) 
that  act  between  suspended  solutes  and  between  solutes  and  substrates,  in  the  LB  method. 
Furthermore,  the  governing  equations  used  in  the  LB  method  naturally  take  into  account  both 
viscous  and  inertial  fluid  effects;  consequently,  the  LB  method  can  be  used  to  study  finite 
Reynolds  number  flows.  Additionally,  the  viscous  term  can  be  modified  to  take  into  account 
non-Newtonian  fluid  effects.  To  date,  microfluidic  simulation  capabilities  deal  only  with  the 
suspending  fluid.  The  critical  need  is  for  a  simulation  capability  that  can  accurately  predict  the 
dynamic  behavior  of  the  particles  in  particle-laden  flows  while  properly  accounting  for  all  of  the 
relevant  physics  encountered  in  microfluidic  systems. 

During  the  program,  we  developed  an  enhanced,  modular  lattice  Boltzmann  (LB) 
simulation  capability  to  predict  the  dynamic  transport  properties  of  particle-laden  flows  in 
microfluidic  sub-systems.  The  specific  physics  capability  that  we  developed  included 
dielectrophoresis  (DEP),  buoyancy,  double  layer  effects  and  electrokinetic  phenomena.  These 
foci  where  driven  by  customer  needs  within  the  BioFlips  program.  Through  the  program  we 
developed  dynamic  simulation  capabilities  to  study  and  characterize  DEP  forces,  DEP  trapping. 
Traveling  wave,  twDEP,  DEP  force  estimation  capabilities,  new  methodologies  to  predict  double 
layer  effects  on  DEP  forces,  and  fluid  and  suspension  effects  in  microchannels  and  sudden 
contractions.  The  capabilities  developed  enable  the  characterization  of  multi-species  DEP 
separations,  and  multi-species  suspension  behavior  in  microenvironments.  The  capabilities  are 
3D  and  parallel.  Additionally,  the  particles  or  suspended  species  in  solution  are  fully  coupled, 
that  is  to  say,  the  particles  interact  with  each  other  and  with  their  environment  without 
approximation.  In  contrast,  most  other  methods,  including  other  methods  supported  within  the 
program,  have  to  use  an  approximation  to  take  into  account  particle-particle  and  particle- 
environment  interactions. 

As  mentioned,  our  effort  was  focused  on  the  overlap  between  the  needs  of  our  customers 
within  the  program  and  our  proposed  goals.  Collaborations  were  formed  with  predominantly  two 
Bioflips  teams,  that  is,  MD  Anderson  Cancer  Center  (Dr.  Peter  Gascoyne)  and  UC  Davis  ( 
Professor  Rosemary  Smith).  In  addition  to  these  collaborations,  we  also  performed  collaborative 
research  with  Professor  Dorian  Eiepmann  (UC  Berkeley)  and  had  some  interactions  with 
Professors  Carl  Meinhart  and  Juan  Santiago.  Through  the  program,  we  are  still  supporting 
Professor  Rosemary  Smith’s  BioFlips  efforts  today,  and  we  have  now  co-authored  two  new 
proposals  with  Professor  Gascoyne.  Also  as  a  result  of  the  program,  we  have  hired  one  of 
Professor  Eiepmann’ s  former  students.  Additionally,  Professor  Meinhart  and  I  have  met  several 
times  to  discuss  research,  and  I  have  now  sponsored  a  student  of  his  to  work  with  me  at  EENE 
twice.  I  have  also  visited  Professor  Santiago  and  given  an  invited  talk  at  Stanford.  I  continue  to 
support  and  mentor  a  student  of  his  that  works  at  EENE. 


The  lattice  Boltzmann  capability  is  applicable  to  a  wide  variety  of  phenomena  relevant  to 
microfluidics  and  other  fields  that  involve  complex  interactions  in  particle-laden  flows.  The 
proposed  capability  is  revolutionary  in  that  it  is  focused  on  the  individual  and  collective  behavior 
of  macromolecules  for  particle-laden  flow  in  microfluidic  devices.  Because  the  relevant  physics 
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involved  will  be  accounted  for,  the  proposed  capability  will  be  useful  in  revealing  new  physical 
phenomena  that  could  in  turn  lead  to  new  microfluidic  technologies.  Furthermore,  this  capability 
constitutes  a  general  framework  where  the  user  can  include  field  and  force  effects  of  choice,  e.g., 
trapping  of  magnetic  beads  via  magnetic  fields,  acoustic  focusing  of  macromolecules,  electro- 
osmotic  effects,  electrophoretic  and  magnetophoretic  separations  to  name  a  few.  Beyond  the 
microfluidics  community,  it  has  wide  application  and  will  enable  the  study  of 
colloidal/macromolecular  transport  in  physiological  systems,  such  as,  blood  filtration  in  the 
kidney,  contaminant  transport  in  soils,  and  it  will  contribute  to  understanding  new  physics 
associated  with  nanoparticle  behavior  in  support  of  the  Nanotechnology  initiative.  This 
revolutionary  dynamic  simulation  capability  will  provide  device  designers  a  valuable  tool  to 
augment  the  device  design  process  and  to  explore  new  sub-system  design  concepts  via  computer 
simulation.  The  proposed  simulation  capability  will  include  the  necessary  physics  modules  to 
accurately  account  for  external  field  effects,  intermolecular  force  interactions,  and  phoretic 
transport  effects  on  individual  and  collections  of  macromolecules  in  microchannels  (particle 
laden  flows). 

Methods,  Assumptions,  Procedures 
Lattice  Boltzmann  Method 

The  LB  method  is  a  mesoscopic  description  of  fluid  behavior.  In  this  approach,  one  does  not 
track  individual  fluid  molecules  but  tracks  a  probability  distribution  function  that  represents  the 
collective  behavior  of  fluid  molecules  in  local  regions.  The  governing  equation  solved  in  the 
lattice  Boltzmann  method,  i.e.,  the  lattice  Boltzmann  equation  (LBE),  is  a  discrete  form  of  the 
Boltzmann  transport  equation  [1].  The  following  equation  is  the  LBE  in  the  simplest  form: 

f.{x  +  ^x,t  +  ^t)=fXx,t)+  (31) 

T 


fi  is  the  z*  component  of  the  single -particle-velocity-distribution  function,  r  is  the  scalar 
relaxation  time,  which  is  directly  related  to  the  kinematics  viscosity,  v,  of  the  fluid,  i.e., 
v'=(2r-l)/6,  and is  the  equilibrium  distribution  function.  The  functionality  of  the 

equilibrium  distribution,  f‘’‘‘{x,t),  in  the  last  term  of  Eq.  (3.1)  is  chosen  to  ensure  that  both  mass 

and  momentum  are  conserved.  Eurthermore,  fi‘"^{x,t)  is  quadratic  in  the  fluid  velocity  and 
therefore  accurately  captures  fluid  inertial  effects.  In  the  limit  of  small  Mach  and  Knudsen 
numbers,  one  can  perform  a  Taylor  series  like  expansion  on  Eq  (3.1).  and  recover  the  Navier- 
Stokes  equations  of  motion;  therefore,  this  approach  is  directly  applicable  to  microfluidic 
systems. 

The  single  particle  distribution  function  is  a  function  of  “lattice”  velocity,  z,  position,  x,  and  time, 
t.  Eq.  (3.1)  is  evolved  on  the  lattice.  /  is  a  moment  baring  function  and  contains  all  of  the 
important  fluid  information  locally  at  each  lattice  site.  The  zero*  discrete  moment. 
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p 


(3.2) 


gives  the  local  fluid  density,  where  p  is  the  fluid  density.  The  first  discrete  moment  with  respect 
to  the  lattice  velocity,  eu  gives  the  local  momentum  density, 

(3-3) 

i 

Finally,  the  velocity  distribution  function  can  be  used  to  calculate  the  force  and  torque  acting  on 
a  stationary  or  mobile  objects  (of  desired  geometry)  using  Ladd’s  [2,  3]  half-link  model. 

To  enforce  the  no-slip  condition  at  solid  surfaces,  the  LB  approach  uses  the  simple  bounce-back 
scheme.  Specifically,  if  the  component  of/ is  advected  from  a  given  fluid  lattice  site  toward  a 
solid  phase  lattice  site  then  that  component  is  reflected  back  to  where  it  originated  and  collides 
with  other  incoming  components  of  the  distribution  function  at  the  originating  lattice  site.  In  this 
way,  the  presence  of  the  solid  phase  is  propagated  back  to  the  bulk  fluid. 

The  bounce-back  condition  affords  an  extreme  amount  of  flexibility.  Specifically,  the  solid  phase 
is  accounted  for  in  a  logical  array  that  includes  all  lattice  sites  in  the  system.  If  a  lattice  site  is 
“true”,  i.e.,  solid  phase,  the  bounce  back  condition  is  enforced.  Otherwise,  it  is  not  enforced.  As 
a  result,  it  is  very  easy  to  include  walls,  obstacles  and  mobile  particles  into  the  capability.  The 
user  is  only  constrained  by  their  ability  to  mathematically  describe  the  object,  i.e.,  to  flag  the 
appropriate  lattice  sites. 

The  lattice  Boltzmann  capability  was  constructed  in  a  modular  fashion.  This  permitted  the 
inclusion  of  proposed  physics  modules.  The  modules  included  in  this  work,  however,  were 
chosen  because  they  represented  main  categories  of  physical  phenomena  that  are  of  primary 
importance  to  biological  microfluidics,  i.e.,  external  field  effects,  intermolecular  forces, 
electrokinetic  phenomena,  and  non-Newtonian  fluid  models.  Shown  in  the  Figure  2  below  is  a 
diagram  of  the  LB  capability  architecture. 


While  the  LB  method  accurately  handles  all  of  the  fluid  mechanics,  it  calculates  everything  in 
“lattice  space.”  As  a  result,  hydrodynamic  forces  in  lattice  space  need  to  be  made  dimensionless 
with  the  appropriate  lattice  parameters  and  converted  to  real  space  force  units.  This  type  of  force 
conversion  is  commonplace  in  the  lattice  Boltzmann  community.  In  the  proposed  work  here,  we 
will  develop  a  set  of  conversion  factors  to  convert  real  space  forces  to  lattice  space  forces: 

(3.4) 

Here  F  lb  is  the  force  in  lattice  space,  is  the  force  conversion  kernel  and  Fy  is  the  force 
from  the  module  of  interest  in  real  space.  This  will  permit  the  summing  of  all  force  contributions 
experienced  by  each  particle  in  the  flow  configuration.  During  the  initial  stages  of  this  work,  we 
verified  this  simple  relationship  on  a  well  known,  controllable  problem,  namely,  for  a  sphere 
sedimenting  in  a  quiescent  fluid  (buoyancy  force). 
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Figure  2  A  basic  block  diagram  showing  tbe  architecture  of  the  module  LB 
capability.  The  force  filter,  where  time  and  length  scales  are  converted 
from  real  space  to  lattice  space  and  from  lattice  space  to  real  space,  is 
central  to  the  overall  capability. 

External  fields/Body  forces/Dielectrophoresis 

In  the  presence  of  non-uniform  electric  fields,  macromolecules  are  either  attracted  to  local 
maxima  or  local  minima  in  the  E-field  strength.  This  can  be  useful  for  separating  biological 
macromolecules  by  collecting  them,  i.e.,  positive  DEP,  or  separating  them  by  physical  properties 
through  levitating  macromolecules,  i.e.,  negative  DEP,  into  a  fluid  flow  field.  Eor  parallel 
electrodes  orthogonal,  see  Eigure  3  below,  to  the  flow  direction  placed  at  the  base  of  a 
rectangular  channel,  the  macromolecules  for  positive  DEP  are  be  attracted  to  electrode  strip 
edges  where  the  E-field  strength  is  a  maximum. 


Figure  3  Left  panel:  Interdigitated  gold  electrodes.  The  flow  is 
from  left  to  right.  The  e-field  strength  is  at  a  maximum  at  the 
electrode  edges.  The  right  panel  shows  the  e-field  strength  at  the 
edge  of  the  electrodes  as  a  function  of  applied  voltage. 
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For  negative  DEP,  particles  are  repelled  (levitated)  from  regions  of  high  field  strength  up  into  the 
flow  field.  Depending  on  the  particle  properties,  differing  species  will  be  levitated  to  different 
positions  in  the  flow  field  and  therefore  will  separate  from  each  other.  This  is  known  as  field 
flow  fractionation,  or,  FEE.  The  simplest  form  of  the  DEP  force  is  given  by  the  following 
expression: 


E-dep  )  -  '^■ns^a  Re 


s  -s 
2 _ 1_ 

e  +2^ 

V  2  1  y 


Y.{e-e). 


(4.1) 


Where  a  is  the  particle  radius,  s*  is  the  complex  permittivity  of  the  phase  andE®  is  the  applied 

electric  field  intensity.  The  term  in  parenthesis  is  known  as  the  Clausius-Mossotti  factor.  It  is 
this  factor,  based  on  fluid  properties  and  the  frequency  of  the  electric  field  that  determines 
whether  one  has  positive  or  negative  DEP.  This  force  is  calculated  in  real  space  and  passed 
through  the  force  filter  and  converted  to  lattice  space  to  be  summed  with  the  forces  calculated 
using  the  EB  method.  While  the  DEP  force  is  small  in  the  “macroscopic”  flow  world,  it  is  very 
important  and  useful  in  the  microscope  world  of  Bio-flips.  Furthermore,  The  DEP  flow  cell  to 
be  studied  in  this  proposal  represents  a  fully  integrated  microfluidic  subsystem  that  is  relevant  to 
several  researchers  and  is  proposed  here  to  demonstrate  the  EB  capability  on  this  fully  integrated 
microsystem. 

To  model  such  micro-unit  operation  we  define  a  global  coordinate  and  local  coordinate  systems, 
see  Figure  4  below. 


Periodic  Simulation  GeU 


Figure  4  Typical  Dielectrophoretic  separator.  The  dashed  rectangular  region  represents  the 
simulation  cell.  The  simulation  cell  moves  with  the  target  species  in  the  flow  direction  toward  the 
interdigitated  electrode  array,  denoted  hy  gold  rectangles. 

Here  the  DEP  separator  is  described  in  full  in  the  global  coordinate  system  and  the  target  species 
is  tracked  in  the  local  simulation  cell  coordinate  system,  the  simulation  cell  is  denoted  in  dashed 
lines  with  the  particle  at  the  centroid.  The  simulation  cell  translates  with  the  particle  as  it 
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approaches  the  electrode  array.  The  simulation  cell  has  periodic  boundary  conditions;  therefore, 
this  “single”  particle  simulation  set-up  actually  takes  into  account  dilute  concentration  effects 
through  the  target  species  periodic  images. 


Intermolecular/Colloidal  force  interactions 


Intermolecular/Colloidal  force  interactions  act  in  a  pair  wise  fashion  (two  body  near-field 
interactions);  therefore,  when  particles  are  within  sufficient  range  of  each  other  (<  1  mm  for 
electrostatics  and  <100  nm  for  van  der  Waals)  we  sum  their  force  contributions  with  the  forces 
calculated  by  the  LB.  Specifically,  we  use  DLVO  theory  to  calculate  these  forces  pass  them 
through  the  force  filter,  see  Eq.  (3.4),  and  sum  the  result  with  the  calculated  LB  force. 


The  approach  to  incorporate  these  forces  is  precisely  the  same  approach  that  the  PI  uses  to 
correct  for  near-field  hydrodynamic  interactions  for  lubrication  forces.  Following  the  same 
approach,  the  non-hydrodynamic  forces  are  summed  with  the  Pi’s  hydrodynamic  force 
correction  to  get  the  total  force  correction  for  a  particular  particle: 


~SEf 
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(5.1) 


where  the  superscripts  T,  H,  and  NH  represent  the  total  force,  the  near-field  hydrodynamic  force, 
and  non-hydrodynamic  force  corrections  respectively.  This  generalized  framework  permits  the 
rapid  inclusion  of  additional  intermolecular/colloidal  forces. 


All  biological  macromolecules  experience  intermolecular  force  interactions  with  each  other  in 
solution,  e.g.,  repulsive  electrostatic  and  attractive  van  der  Waals  force  interactions.  These  force 
interactions  can  cause  non-intuitive  particle  behavior  particularly  when  there  are  many 
interacting  particles  in  confined  domains,  e.g.,  particle  aggregation,  wall  adhesion,  etc.  In  micro¬ 
scale  channels  these  force  interactions  become  increasingly  more  important.  Electrostatic 
interactions  act  on  the  length  scale  of  the  Debye  screening  length,  which  is  a  measure  of  the 
thickness  of  the  electric  double  layer,  i.e.,  electrostatic  force  interactions  act  over  a  distance 
usually  less  than  one  micron.  Where  the  Debye  screening  length,  is  inversely  proportional  to 
the  square  root  of  the  ionic  strength  of  the  solution: 

=  cl  ^  Electrolyte  Concentration  ,  (5.2) 


where  C  is  a  constant  dependent  on  the  salt  species  in  solution  [4].  The  units  on  Eq.  (5.2)  are 
nanometers,  nm.  The  attractive  van  der  Waals  forces  on  the  other  hand  act  on  the  length  scale  of 
100  nm  or  less  [4].  Because  both  of  these  force  interactions  act  in  a  pair  wise  fashion  (two  body 
near-field  interactions),  their  presence  is  accounted  for  in  form  of  an  interaction  potential,  which 
provides  a  framework  for  inclusion  of  other  pair  wise  force  interactions  that  are  of  interest. 
These  force  interactions  are  present,  and  in  micro-flows,  they  must  be  taken  into  account. 
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Phoretic  transport  mechanisms 


Electrophoresis  and  electro-osmosis  transport  mechanisms  are  very  common  in  all  microfluidic 
systems.  The  mechanism  that  causes  each  phenomenon  is  the  same  in  both  cases.  Specifically,  a 
charged  surface  in  an  electrolyte  solution  has  a  cloud  of  counter  ions  that  seeks  to  electrically 
neutralize  the  surface,  or,  achieve  electro-neutrality.  The  charged  surface  and  the  cloud  of 
counter  ions  are  what  is  known  as  the  electric  double  layer.  When  subject  to  a  uniform  field, 
usually  a  DC  field,  the  electrode  with  the  same  charge  as  the  counter  ion  cloud  causes  a  flow  by 
repelling  the  like  ions.  For  stationary  surfaces  this  gives  rise  to  electro-osmotic  flows,  and  for 
mobile  surfaces,  this  gives  rise  to  electrophoretic  motion  of  the  mobile  species.  The  velocity  of 
counter  ion  repulsion  can  be  described  as  a  slip  velocity  at  the  charged  surface  [5].  For 
electrophoresis  and  electro-osmosis  the  slip  velocity  is  given  by 

(6.1) 


Where  Ms  is  the  slip  velocity  at  the  surface,  s  is  the  permittivity  of  the  suspending  fluid,  C,  is  the 
zeta  potential  (the  measured  surface  potential  with  adsorbed  counter  ions),  and  is  the  applied 
electric  field.  This  velocity  is  the  electrophoretic  velocity  at  which  a  particle  with  zeta  potential 
C,  travels  through  a  fluid  of  permittivity  of  s  and  viscosity  of  //  subject  to  a  uniform  electric  field 
evaluated  at  the  surface,  E^.  For  a  stationary  surface,  Eq.  (6.1)  describes  the  boundary  condition 
at  a  charged  surface  that  drives  electro-osmotic  flow.  Other  phoretic  transport  phenomena,  such 
as,  diffusiophoresis  and  thermophoresis,  can  also  be  described  in  terms  of  a  slip  velocity  [5]; 
therefore,  this  represents  a  general  framework  for  including  this  class  of  transport  phenomena. 

While  the  physics  above  appears  somewhat  disjoint,  their  sum  represents  the  coupled  phenomena 
that  exist  in  many  real  Bio-fluidic-chips.  For  example,  it  is  not  uncommon  to  use 
dielectrophoresis  to  trap  or  sort  biological  macromolecules.  In  these  micro-systems,  bio¬ 
molecules  convect  with  the  fluid  in  microchannels.  The  bio-molecules  in  particle-laden  flow 
have  three  predominant  effects:  External  fields,  i.e.,  gravity  and  AC  electric  fields,  are 
incorporated  into  the  FB  method  through  analytic  and  semi-analytic  expressions  [6].  The  forces 
induced  by  these  two  external  field  effects,  buoyancy  and  DEP  forces,  will  also  be  passed 
through  the  force  filter  and  summed  with  the  force  as  predicted  by  the  FB  method,  see  Eq.  (3.4) 
above.  Additionally,  the  resulting  capability  is  3D  and  parallel,  using  F90  and  MPI. 

Results  and  Discussion 

The  resulting  effort  was  directed  by  the  needs  of  the  BioFlips  Pis;  therefore  we  focused  in  the 
three  following  areas: 

1 .  Dielectrophoretic  separations  achieved  by  preferential  particle  trapping, 

2.  Particle  transport  independent  of  fluid  motion,  traveling  wave  Dielectrophoresis,  and 

3.  Suspension  properties  in  micro  flows,  sudden  contraction  in  micro  flows. 
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In  addition  to  these  three  thrust  areas,  we  explored  some  anomalous  effects  experience  in  DEP 
separations.  In  particular,  we  explored  and  elucidated  the  effects  of  the  electric  double  layer  on 
the  polarizability  of  the  target  species. 

Dielectrophoresis 

As  discussed  above.  Dielectrophoresis,  DEP,  arises  as  a  result  of  an  induced  dipole  in  a  non- 
uniform  electric  field,  which  results  in  a  net  Eorenz  force  either  levitating  or  attracting  the  target 
species  [7].  Several  Pi’s  in  the  BioElips  program  employ  DEP  to  perform  micro-separations  and 
sample  preparation.  Shown  in  Eigure  5  below  is  the  time  average  DEP  force  for  standing  wave 
DEP. 


Figure  5  Interdigitated  electrode  configuration  and  the  electric  field  intensity  produced  by  this  electrode 
configuration.  The  gold  rectangles  in  the  right  panel  represent  electrodes  in  the  array.  The  surface  plot  of 
electric  field  intensity  emanating  from  the  electrode  surface. 

The  DEP  force  is  directly  proportional  to  the  gradient  in  the  electric  field  intensity,  ECE  ,[8]  as 
shown  in  the  right  panel  of  Eigure  5.  Given  that  the  electric  field  intensity  is  a  maximum  at 
electrode  edges,  the  DEP  force  is  also  a  maximum  at  electrode  edges,  see  Figure  6. 

The  gold  rectangles  represent  the  electrodes  and  the  z  =  0  surface  or  floor  of  the  electrode  array. 
Note  that  the  gradients  are  indeed  much  stronger  at  the  electrode  edges.  The  z-direction  or  X3- 
direction  is  always  negative,  the  vector  pointing  to  the  electrode  edge.  When  the  Clausius- 
Mossotti  factor  is  “positive”  we  have  positive  DEP  and  the  target  species  is  attracted  to  regions 
of  high  field  intensity,  toward  the  electrode  edge.  When  the  Clausius-Mossotti  factor  is 
“negative”  we  have  negative  DEP  and  the  target  species  is  repelled  or  levitated  due  to  the 
presence  of  the  applied  field. 
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To  assist  designers,  we  explored  the  effects  of  changing  the  electrode  width  and  spacing  on  the 
magnitude  of  the  gradient  in  the  electric  field  intensity,  [8],  see  Figure  7 
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Figure  6  Gradients  in  the  electric  field  intensity.  The  left  panel  is  the  gradient  in 
the  flow  direction  and  the  right  panel  is  the  gradient  in  the  vertical  direction. 
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Figure  7  The  dependence  of  the  magnitude  of  the  gradient  in  the  electric  field  strength.  The 
left  panel  is  the  dependence  on  the  applied  voltage.  The  right  panel  is  the  dependence  on  the 
spacing  between  electrodes. 
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As  shown  the  gradient  in  the  eleetrie  field  intensity  steepens  as  the  eleetrode  width  is  reduced. 
Also  shown  on  the  right  hand  panel,  the  gradient  in  the  electric  field  intensity  increases  as  the 
electrodes  are  brought  closer  together.  For  further  details  the  reader  is  referred  to  [8]. 

To  determine  the  critical  Reynolds  number  of  operation  for  a  given  applied  voltage,  we 
performed  a  serious  of  numerical  experiments  to  construct  a  phase  diagram,  see  Figure  8  below. 
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Figure  8  Determination  of  the  critical  Reynolds  number  for  a  given  electrode  configuration,  applied 
voltage  and  particle  radius.  The  left  panel  shows  particle  translation  as  a  function  of  Reynolds  number 
relative  to  an  electrode,  denoted  in  gold.  The  right  panel  is  a  phase  diagram  for  predicting  particle 
trapping  given  a  Clausius-Mossotti  factor  and  a  channel  Reynolds  number. 


To  assist  designers  to  know  what  flow  rates  or  Channel  Reynolds  numbers  that  they  can  operate 
their  separators,  we  performed  a  series  of  simulations  to  determine  the  critical,  channel  Reynolds 
number  as  a  function  of  Clausius-Mossotti  factor  for  target  species  trapping.  By  doing  a  series  of 
experiments,  we  were  able  to  construct  a  phase  diagram,  see  right  panel,  for  rapid  estimation  of 
DEP  separator  operating  conditions,  [9,  10]. 
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Figure  9  Critical  Reynolds  number  for  particle  trapping  as  a  function  of  applied  voltage  (left  panel)  and  as  a 
function  of  target  species  radius,  a  (right  panel). 


In  addition  to  these  studies,  we  also  validated  our  dynamie  simulation  eapability  with  theory  by 
predicting  the  dependence  of  the  critical  Reynolds  number  on  the  square  of  the  applied  voltage 
and  the  target  species  radius.  Our  simulations  exhibited  the  expected  dependence.  Two  additional 
parameters  that  practitioners  need  to  know  are  1)  how  long  does  the  electrode  array  need  to  be  to 
and  2)  what  depth  of  channel  should  be  used  for  an  applied  voltage  to  ensure  target  species 
manipulation/trapping;  therefore,  we  performed  a  series  of  numerical  experiments  to  study  the 
depth  of  penetration  of  the  Dielectrophoretic  force  in  a  DEP  separator,  see  Figure  10. 
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Figure  10  Dielectrophoretic  particle  trapping  as  a  function  of  particle  height  for  a  given 
applied  voltage,  electrode  configuration,  Clausius-Mossoti  factor  and  species  radius. 
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For  the  conditions  show  on  the  right  side  of  Figure  7  the  DEP  field  reaches  approximately  40  um 
into  the  fluid.  Therefore  in  the  design  of  a  separator  that  operates  at  similar  operating  conditions, 
the  user  will  need  to  adjust  the  depth  of  the  separator  channel  and  the  length  of  the  overall 
electrode  array  to  ensure  that  the  species  of  interest  is  trapped  with  high  efficiency  out  of  the 
sample. 

To  make  this  unique  capability  more  useful,  we  extended  it  to  include  multiple  species.  This  was 
as  a  result  of  our  collaboration  with  MD  Anderson  cancer  center.  In  our  support  of  their  effort, 
we  developed  capabilities  to  predict  species  cross  over  frequencies  and  to  predict  suspension 
effects  on  separator  performance.  In  the  MD  Anderson  application  they  were  particularly 
interested  in  separating  the  various  white  cells  from  whole  blood  samples  to  monitor  changes  in 
concentration  in  response  to  exposure  to  pathogen.  The  cross  over  frequencies,  where  the 
Clauisus-Mossotti  factor  is  zero,  for  T  and  B  Lymphocytes,  Monocytes  and  Granulocytes  are 
shown  in  Figure  1 1 . 


Figure  11  Dielectrophoretic  response  of  white  blood  cells. 

The  goal  is  to  determine  the  cross  over  frequencies  to 
optimize  the  Dielectrophoretic  separation  protocol. 

By  knowing  the  cross  over  frequencies,  the  DEP  separator  protocol  can  be  optimized  to  effect 
desired  separations.  Finally,  we  developed  a  multi-species  DEP  simulation  capability,  see  Figure 
12. 

In  this  unique  capability,  which  uses  the  simulation  configuration  as  described  in  Figure  4, 
researchers  can  now  explore  separation  protocols  as  a  function  of  the  species  present,  fluid 
conditions  and  applied  voltages. 

With  this  capability,  we  explored  the  design  of  a  traveling  wave  DEP  separator  with  Professor 
Rosemary  Smith  at  UC  Davis.  Shown  below  in  Eigure  13  is  a  cartoon  of  their  separation 
chamber. 
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•  Multiple  species  with  unique  properties: 

-  radius, 

-  density, 

-  permittivity, 

-  condnctivity, 

-  fluid  permittivity, 

-  fluid  conductivity, 

-  variabie  field  strength  and  freqnency, 

-  electrode  width  and  spacing, 

-  buoyancy  effects. 


Figure  12  Multispecies  simulatiou  capability.  Species  iu  solutiou  cau  be  attributed 
with  uuique  physical  properties,  e.g.,  Dielectrophoretic  properties. 
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Figure  13  Traveling  wave  Dielectrophoretic  separator. 
Antibody  analyte  aggregates  translating  independent  of 
the  fluid  to  a  detection  region. 


In  the  above  figure,  a  sample  of  body  fluid  is  drawn  into  a  chamber  of  quiescent  fluid  and 
liposomes  form  aggregates  mediated  by  target  analyte  in  solution.  The  liposome  aggregates  are 
translated  independent  of  fluid  flow  by  traveling  wave  DEP  to  a  detection  region.  The  traveling 
wave  force  is  dependent  on  both  the  standing  wave  and  the  gradient  in  the  phase  of  the  applied 
field  as  shown  in  Figure  14. 
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Figure  14  Traveling  wave  Dielectrophoresis  translate  the  target  species 
independent  of  the  fluid  (top).  The  driving  force  is  the  gradient  in  the  phase 
of  the  electric  field  (right  panel). 

Note  that  the  time  averaged  force  contains  both  the  standing  and  traveling  wave  components; 
therefore,  the  target  species  is  levitated  and  translated  during  traveling  wave  DEP.  Also  note  that 
the  traveling  wave  piece  is  dependent  on  the  imaginary  part  of  the  Claus ius-Mossotti  factor, /cm- 
To  calculate  the  translocation  velocity,  we  determined  the  particle  velocity  that  is  equivalent  with 
the  tw-DEP  force.  In  Eigure  15  below,  a  plot  of  the  translocation  velocity  as  a  function  of  height 
above  the  electrode  array  is  given. 
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Figure  15  Translocation  velocity  of  target  species  predicted  from  the  balance  of  the  traveling 
wave  force  and  the  viscous  force  exerted  hy  the  fluid  on  the  body.  The  right  panel  is  a  map  of 
the  resulting  traveling  wave  velocities  as  a  function  of  position  from  the  electrode  array. 
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As  shown  by  Green  et  al.,[l  1]  there  exists  an  optimal  height  above  the  electrode  for  traveling 
wave  DEP.  This  critical  height  is  when  the  ratio  of  the  height  with  the  electrode  width  is  ~1.  In 
direct  support  of  the  UC  Davis  effort,  we  determined  translocation  velocities  for  various  particle 
sizes  and  applied  voltages,  see  Figure  16  below. 
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Figure  16  Predicted  translocation  velocities  as  a  function  of  particle  size  and  applied 
voltage. 

In  this  study  we  chose  a  physiological  solution  conductivity  and  calculated  translocation 
velocities  at  various  aggregate  radii  for  two  applied  voltages.  These  studies  directly  impacted  the 
UC  Davis  device  design. 

Double  layer  effects 

A  common  phenomenon  kept  reappearing  in  our  support  of  both  MD  Anderson  and  UC  Davis 
(as  well  as  other  efforts  that  use  DEP).  This  phenomenon  was  that  target  species  tended  to 
exhibit  the  opposite  DEP  behavior.  More  specifically,  the  calculated  Clausius-Mossotti  factor 
would  predict  the  sign  of  the  DEP  force  acting  on  the  target  species,  but  the  species  would 
respond  as  if  the  Clausius-Mossotti  factor  had  an  opposite  sign  (while  the  magnitude  of  the  force 
was  not  necessarily  equal  and  opposite  to  what  was  predicted;  therefore,  we  looked  into  this 
phenomenon.  We  noted  that  while  DEP  is  a  volumetric  force  that  the  Clausius-Mossotti  factor 
depended  heavily  on  the  surface  properties  of  the  target  species;  furthermore,  we  observed  that 
the  species  that  exhibited  this  phenomenon  had  a  surface  charge  and  therefore  an  electric  double 
layer.  This  led  us  to  develop  a  new  approach  for  calculating  the  Clausius-Mossotti  factor  to  take 
into  account  not  only  the  Stern  layer  but  the  diffuse  double  layer  as  shown  in  Figure  17. 

The  electric  double  layer  consists  of  two  parts,  1)  the  bound  counter  ions  on  the  surface,  the 
Stern  layer  and  2)  the  diffuse  counter  ion  cloud  that  seeks  to  neutralize  the  charge  of  the  surface, 
a  condition  known  as  electro-neutrality.  By  calculating  the  conductivity  and  the  change  in 
permeability  in  the  double  layer  and  using  the  shell  model  [7],  we  incorporate  the  presence  of  the 
electric  double  layer  when  predicting  the  effective  properties  of  the  target  species.  As  it  turns 
out,  all  biological  species,  e.g.,  analytes,  and  polystyrene,  which  is  ubiquitous  in  microfiuidics  all 
have  a  surface  charge  and  therefore  exhibit  anomalous  response  to  DEP  forces.  Therefore  such 
an  understanding  and  predictive  capability  is  critical  to  the  successful  implementation  of  DEP  to 
effect  desired  separations.  Many  of  those  who  are  new  to  the  use  of  DEP  erroneously  just  switch 
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the  sign  on  the  Clauisius-Mossotti  factor  to  account  for  the  unexpected  change  in  behavior.  This 
has  occurred  even  within  the  BioFlips  efforts  when  predicting  the  DEP  response  of  polystyrene 
beads. 


Figure  17  The  details  of  the  electric  double  layer  iucorporated  iuto  a 
respeseutative  particle  with  the  appropriate  effective  couductivity,  a,  aud 
effective  dielectric,  e,  properties. 


In  Figure  18  below,  typical  properties  for  polystyrene  beads  used  in  biological  applications  are 
given.  As  shown,  we  assumed  a  radius  of  216  nm.  Also  shown  is  a  cartoon  of  the  electric  double 
layer.  Given  these  properties,  we  applied  our  new  capability  to  predict  changes  in  the  Clauisus- 


Radius  =216  nm 
cTp  ~  10-^2  ms/m 

=  2.55 

t;  ~  90  mv 

Hughes  et  al.,  1999 


Diffuse  layer 


Figure  18  Polystyreue  head.  Physical  properties  takeu  from  Hughes  et  al.,  1999. 

Mossotti  factor  as  a  function  of  ion  strength  of  the  solution.  The  results  for  three  solution 
conductivities,  concentrations  of  [NaCl]  are  given  below  in  Figure  19. 
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Constant  Zeta  potential,  variation  in  concentration 


Figure  19  The  real  part  of  the  frequency  dependent  Clausius- 
Mossotti  factor  for  Polystyrene  as  a  function  of  ionic  strength 
of  the  solution. 


Without  accounting  for  the  electric  double  layer,  one  would  always  predict  a  negative  Clausius- 
Mossotti  factor  and  therefore  would  always  predict  a  repulsive  DEP  force.  In  practice  however, 
polystyrene  beads  are  “weakly”  trapped  depending  on  the  ionic  strength  of  the  solution.  As 
shown  above  in  Figure  19,  our  new  approach  predicts  the  “positive”  DEP  as  a  result  of  double 
layer  effects.  Also  note  that  the  “cross-over”  frequeney  deereases  with  inereasing  solution 
eonduetivity. 

To  explore  the  effects  on  simulants,  such  as,  Avidin,  we  also  made  predictions  to 
demonstrate  the  new  capability  on  proteins,  see  Figure  20. 

Protein  ^ize  Range 

-  10  nm 


Zeta  Potential 

•  Dependent  on  pH 

•  <^G  [60...-60mv] 

(Wiqcek  and  Chibowski,  2002) 


Protein 


Figure  20  Typical  properties  of  proteins. 
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In  is  important  to  note  that  the  Zeta  potential,  whieh  is  extremely  important  in  estimating  the 
conductivity  and  permittivity  of  the  electric  double  layer,  is  highly  sensitive  to  the  solution  pH. 
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Figure  21  The  real  part  of  the  Clausius-Mossotti  factor  for  a  typical  protein 
at  two  different  fluid  conductivities. 


Again,  in  Figure  21,  we  see  the  dramatic  effects  of  the  electric  double  layer  on  the  DEP  response 
of  a  protein.  Consistent  with  the  typical  prediction  of  the  Clausius-Mossotti  factor,  the  response 
without  taking  into  account  the  effects  of  the  electric  double  layer  always  gives  negative  DEP; 
however,  when  accounting  for  the  electric  double  layer,  the  response  is  significantly  different, 
and  in  the  case  shown,  the  response  now  always  predicts  positive  DEP.  We  preformed  this 
analysis  specifically  for  the  UC  Davis  team  because  they  were  using  Avidin  as  their  simulant. 
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Figure  22  Liposomes  with  a  mauufactured  hiolayer. 

A  liposome  is  simply  a  spherical  volume  of  fluid  enclosed  by  a  bilayer. 
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Radius  =  SOOnm 
Membrane  Properties 
Capacity  =  6.7mF/m2 


Figure  23  The  real  part  of  the  Clausius-Mossotti  factor  for  a  typical  protein  at 
two  different  fluid  conductivities. 


The  Liposomes  have  very  unique  and  interesting  behavior.  Given  that  the  fluid  inside  the 
liposome  is  identical  to  the  exterior  fluid,  one  would  expect  relatively  uninteresting  behavior.  As 
shown,  however,  the  behavior  is  quite  interesting.  Note  how  the  response  switches  from  negative 
to  positive  DEP  in  the  low  frequency  range  when  the  fluid  conductivity  is  less  than  0.1  mS/m. 
This  elucidates  the  fact  that  while  DEP  is  a  body  force  the  species  response  is  predominantly 
governed  by  the  electrical  properties  of  the  membrane. 

With  this  new  capability,  device  designers  and  practitioners  can  more  accurately  predict  the 
response  of  “charged”  species  and  be  able  to  tune  fluid  properties  to  achieve  desired  separations. 

Suspension  Mechanics 

One  of  the  beauties  of  the  lattice  Boltzmann  method  is  the  ability  to  handle  multiple  suspended 
species  with  various  properties.  In  addition,  one  of  the  main  foci  of  the  BioElips  program  had  to 
do  with  sample  collection  and  fluid  delivery.  In  particular,  at  the  South  Carolina  meeting,  there 
was  a  focus  group  on  micro-needles.  Also,  we  had  interacted  with  UCB  on  drug  delivery  and 
microneedles.  It  became  apparent  that  there  were  issues  that  needed  to  be  understood  and  over¬ 
come  to  unleash  the  use  of  microneedles  in  real  world  applications.  In  particular,  in  the  UCB 
case,  they  were  experiencing  microneedle  clogging  at  very  low  volume  fractions  (0.1%),  see 
Eigure  24. 
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Figure  24  Experimental  results  of  particle  coagulation  through  a  sudden 
contraction  (Results  courtesy  of  Professor  Dorian  Liepmann). 


On  the  left  the  length  scales  of  interest  are  shown,  and  on  the  right  is  an  experimental  plot 
showing  the  applied  pressure  as  a  function  of  time  to  elucidate  the  onset  of  clogging  or  particle 
coagulation.  To  assist  the  UCB  team,  we  developed  simulation  capabilities  to  study  the  fluid 
mechanics  and  suspension  behavior  through  a  sudden  contraction.  As  shown  in  Figure  25  below, 
we  developed  a  new  capability  to  apply  a  pressure  gradient  through  the  contraction  to  study  the 
fluid  velocity  fields,  shear-stress  fields  and  pressure  field[12]. 


Figure  25  Lattice  Boltzmauu  represeutatiou  of  suddeu  coutractiou  (3:1  z-directiou,  2:1  y-directiou) 
iu  a  microueedle.  The  left  pauel  represeuts  the  velocity  field  aud  the  right  pauel  is  the  maguitude  of 
the  pressure  field. 
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Due  to  conservation  of  mass  and  volumetric  flow  rate,  we  see  a  significant  increase  in  fluid 
velocity  through  the  contraction.  This  is  not  necessarily  earth  shattering,  but  what  is  important  to 
note  is  that  the  stream-lines  leading  to  the  contraction  exhibit  a  significant  fluid  acceleration  that 
impart  sufficient  momentum  to  the  suspended  species  to  cross  stream  lines.  The  1  micron  beads 
used  had  a  polymer  surface  coating  to  disrupt  the  electric  double  layer,  but  this  surface  coating 
serves  as  a  hydrophobic  corona  that  mediates  particle  coagulation  when  particles  come  close  to 
contact.  As  a  result,  due  to  the  contraction,  the  particles  cross  stream-lines  and  interact  and 
coagulate.  Also  shown  is  the  pressure  field,  it  is  interesting  to  note  that  the  pressure  field  has  a 
much  more  significant  gradient  in  the  narrow  portion  of  the  micro-needle.  While  this  is 
intuitively  obvious,  this  result  graphically  demonstrates  this  physical  phenomenon.  With  the  LB 
method,  we  can  easily  vary  the  chaimel  geometry;  therefore,  this  capability  is  ideal  for  exploring 
and  optimization  of  micro-needle  design. 

In  keeping  with  the  specific  aims  of  the  proposed  effort,  we  developed  the  capability  to 
handle  suspensions  in  microchannels,  from  dilute  to  concentrated  suspensions,  see  Figure  26 
below  for  a  typical  simulation  cell  at  40%  concentration  by  volume. 


Figure  26  Lattice  Boltzmann  capability  to  handle  concentrated  suspensions. 

The  left  panel  is  the  suspension  initialization  at  a  specific  concentration.  The 
right  panel  is  a  snap-shot  of  particle  positions  in  a  rectangular  channel  in 
pressure  driven  flow. 

Also,  shown  is  a  frame  from  a  dynamic  simulation  of  a  suspension  in  a  straight  channel.  Here 
each  particle  position  is  tracked  and  the  particles  are  fully  coupled  with  each  other  and  the 
environment,  that  is  to  say  they  are  interacting  hydrodynamically  with  each  other  and  the 
channel  walls  with  no  approximation.  The  color  in  both  panels  is  insignificant;  it  is  just  an 
artifact  of  MatLab.  Our  capability  not  only  accounts  for  far  field  hydrodynamic  interactions,  but 
we  have  lubrication  interactions  built  in  for  particle-particle  and  particle-wall  interactions. 

To  validate  our  new  capability,  we  compared  our  suspension  results  with  the  gold  standard  put 
forward  by  Nott  and  Brady  [13].  They  discovered  unique  ways  to  characterize  suspension  flow 
between  parallel  walls.  In  keeping  with  their  developments  in  Figure  27  below  we  characterize 
our  results  for  suspension  behavior. 
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Velocity  Profile  Temperature  Profile 


Figure  27  Suspension  velocity  profile  compared  with  a  pure  fluid  velocity  profile  (left  panel),  and  the 
collision  frequency,  e.g.,  suspension  temperature,  on  the  right. 


In  the  left  panel  we  compare  the  suspension  velocity  with  the  velocity  profile  for  pure  fluid.  The 
observed  behavior  is  best  understood  by  taking  into  account  the  finite  size  of  the  suspended 
species.  Note  the  violation  of  the  no  slip  condition  on  the  walls,  x  =  0  and  x  =110.  This  is  due  to 
the  fact  that  the  particle  centroid  cannot  sample  the  no-slip  condition.  Also  note  the  blunting  of 
the  suspension  velocity  profile  across  the  channel.  This  is  also  due  to  the  finite  size  of  the 
suspended  species.  Specifically,  to  get  the  velocity  of  an  individual,  non-interacting  particle,  one 
would  have  to  integrate  the  impinging  velocity  profile  over  the  surface  of  the  entire  particle,  and 
we  have  parabolic  flow,  so  the  velocity  profile  is  non-uniform.  On  the  right  panel,  we  show  an 
interesting  result/metric  used  to  characterize  suspension  behavior  in  confined  domains.  More 
specifically,  we  refer  to  the  suspension  temperature.  This  is  a  temperature  in  the  most 
fundamental  sense,  that  is,  we  characterize  the  numbers  of  particle  collisions  with  other  surfaces 
across  the  channel.  The  majority  of  collisions  occur  where  the  fluid  velocity  gradient  is  steepest. 
In  addition  to  exploring  these  suspension  characteristics,  we  calculated  a  suspension  viscosity  to 
compare  with  the  well-known  Krieger  correlation 


(8.1) 


Here  //  is  the  suspension  viscosity,  is  the  fluid  viscosity,  (j)  is  the  concentration  of  particles, 
and  (j)^  is  the  concentration  of  particles  at  maximum  packing.  Our  result  and  Eq.  (8.1)  were  in 
excellent  agreement.  This  is  extremely  useful  to  device  designers,  that  is,  it  is  paramount  to  know 
the  suspension  viscosity  when  deciding  on  pressure  drops  and  volumetric  flow  rates. 
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Electro-osmotic  flow 


Finally,  we  adapted  the  LB  method  to  handle  electro-osmotic  flows.  This  was  accomplished  by 
incorporating  Eq.  (6.1)  in  to  the  LB  formulation  at  charged  surfaces.  The  results  are  shown 
below  in  Figure  28. 


EO  Flow  (2D  Channel) 


Figure  28  The  velocity  profile  produced  by  Electro-Osmotic 
bouudary  couditious  iu  a  rectaugular  chaouel.  The  bouudary 
couditiou  at  the  chaouel  wall  was  set  by  the  Smolucalski  equatioo  for 
thio  double  layers. 


This  simple  example  demonstrates  the  capability  as  applied  to  a  thin  double  layer  configuration 
in  a  straight  channel.  Given  the  flexibility  of  the  LB  method,  this  result  can  easily  be  extended  to 
more  complex  geometries,  variable  Zeta  potential  and  include  suspended  species.  However, 
during  the  Program,  we  were  asked  to  target  our  efforts  specifically  to  what  our  customers 
needed.  Our  customer  base  was  predominately  associated  with  Dielectrophoresis.  Only  in  the  last 
few  months  did  we  identify  a  customer  who  needed  our  capability  to  study  suspension  mechanics 
in  EO  flow;  unfortunately,  we  had  all  of  our  resources  focused  on  DEP  at  the  time.  As  a  result, 
this  work  was  not  performed  by  our  team. 
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Gantt 


Figure  29  Gantt  Chart  for  the  projeet  aims  as  put  forward  in  the 
original  proposal. 


Conclusions 

The  lattice  Boltzmann  method  is  an  extremely  useful  approach  for  studying  particle  dynamics  in 
microflows.  The  approach  readily  permits  the  inclusion  of  new  physics.  The  limitations  include 
the  following:  1)  the  computational  domain  is  limited  to  sub-systems.  As  a  result,  inlet  and  outlet 
boundary  conditions  need  to  be  specified  using  theory  or  another  numerical  approach;  and  2)  the 
approach  is  limited  to  mm  length  scales  and  micro-second  time  scales;  therefore  it  is  an 
improvement  over  MD  simulations,  but  to  study  systems  of  interest,  the  approach  must  be  made 
parallel  and  requires  significant  computational  facilities,  e.g.,  for  large  problems.  In  spite  of  these 
limitations,  the  PI  and  team  found  that  the  LB  method  is  useful,  and  we  are  continuing  to 
develop  this  capability  to  cover  additional  application  space.  Additionally,  all  of  the  physics 
modules  developed  through  this  program  are  stand  alone,  are  of  great  benefit  and  run  on  a  single 
processor  machine. 


25 


Recommendations 


The  lattice  Boltzmann  method  is  extremely  useful,  and  therefore,  I  recommend  continued 
support  and  use  of  this  approach  in  the  characterization  of  micro  flows.  Given  this  foundation,  I 
recommend  that  an  effort  to  integrate  the  LB  method  to  molecular  level  simulation  and 
macroscale  simulation  be  made.  The  resulting  capability  would  be  a  very  powerful 
computational  tool  that  would  link  the  effects  of  molecular  level  phenomena  to  macroscale 
performance. 
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